1 Introduction

A mixed model with random effects was chosen for this multifactor experiment, and analyzed using the limma package in R. This package implements a linear modeling approach and empiracal Bayes statistics. The limma package with the voom method estimates the mean-variance relationship of the log-counts, generates a precision weight for each observation and enters these into the limma empirical Bayes analysis pipeline.

In this model, clade (levels = Clade1, Clade2, Clade3), physiology (levels = Marine, Freshwater), and experimental condition (levels = 15ppt, 0.2ppt) are fixed effects while species (levels = 14) is considered a random effect.

The raw counts file, generated with NumReads from the salmon (version 0.12.0) quantification tool and summarized with the tximport Bioconductor package (version 1.10.1) in R, can be downloaded from an osf repository with this link, then imported into the R framework.

Samples from species with low numbers of replicates were dropped from the raw counts table (F. zebrinus, F. nottii, F. sciadicus). The raw counts table has the following dimensions (genes x samples).

#dim(counts_design)
#[1] 30471   130

# -----------------------
# Format design and counts matrix
# Drop columns with no data
# -----------------------

design <- counts_design[counts_design$Ensembl == 'Empty',]
#design$type <- c("species","native_salinity","clade","group","condition")
drops <- c("X","Ensembl",
           "F_zebrinus_BW_1.quant","F_zebrinus_BW_2.quant",
           "F_zebrinus_FW_1.quant","F_zebrinus_FW_2.quant",
           "F_notti_FW_1.quant","F_notti_FW_2.quant",
           "F_sciadicus_BW_1.quant","F_sciadicus_FW_1.quant","F_sciadicus_FW_2.quant")
transfer_drops <- c("F_sciadicus_transfer_1.quant","F_rathbuni_transfer_1.quant","F_rathbuni_transfer_2.quant","F_rathbuni_transfer_3.quant",
                    "F_grandis_transfer_1.quant","F_grandis_transfer_2.quant","F_grandis_transfer_3.quant",
                    "F_notatus_transfer_1.quant","F_notatus_transfer_2.quant","F_notatus_transfer_3.quant",
                    "F_parvapinis_transfer_1.quant","F_parvapinis_transfer_2.quant",
                    "L_goodei_transfer_1.quant","L_goodei_transfer_2.quant","L_goodei_transfer_3.quant",
                    "F_olivaceous_transfer_1.quant","F_olivaceous_transfer_2.quant",
                    "L_parva_transfer_1.quant","L_parva_transfer_2.quant","L_parva_transfer_3.quant",
                    "F_heteroclitusMDPP_transfer_1.quant","F_heteroclitusMDPP_transfer_2.quant","F_heteroclitusMDPP_transfer_3.quant",
                    "F_similis_transfer_1.quant","F_similis_transfer_2.quant","F_similis_transfer_3.quant",
                    "F_diaphanus_transfer_1.quant","F_diaphanus_transfer_2.quant",
                    "F_chrysotus_transfer_1.quant","F_chrysotus_transfer_2.quant",
                    "A_xenica_transfer_1.quant","A_xenica_transfer_2.quant","A_xenica_transfer_3.quant" ,
                    "F_catanatus_transfer_1.quant","F_catanatus_transfer_2.quant",
                    "F_heteroclitusMDPL_transfer_1.quant","F_heteroclitusMDPL_transfer_2.quant","F_heteroclitusMDPL_transfer_3.quant")
counts<-counts_design[!counts_design$Ensembl == 'Empty',]
rownames(counts)<-counts$Ensembl
design <- design[ , !(names(design) %in% drops)]
counts <- counts[ , !(names(counts) %in% drops)]
design <- design[ , !(names(design) %in% transfer_drops)]
counts <- counts[ , !(names(counts) %in% transfer_drops)]
#dim(design)
#[1]  5 81
dim(counts)
## [1] 30466    81
gene.names<-rownames(counts)
design[] <- lapply( design, factor)

2 Sample Design Matrix

A matrix was generated using the following model:

~0 + physiology*condition*clade

The random effect of species will be taken into account later.

# --------------------
# design cateogories
# --------------------

species<-as.character(unlist(design[1,]))
physiology<-as.character(unlist(design[2,]))
clade<-as.character(unlist(design[3,]))
condition<-as.character(unlist(design[5,]))
condition_physiology<-as.vector(paste(condition,physiology,sep="."))
condition_physiology_clade <- as.vector(paste(condition_physiology,clade,sep="."))
condition_physiology_clade <- as.vector(paste("group",condition_physiology_clade,sep=""))
cols<-colnames(counts)
ExpDesign <- data.frame(row.names=cols,
                        condition=condition,
                        physiology = physiology,
                        clade = clade,
                        species = species,
                        sample=cols)
ExpDesign
design = model.matrix( ~0 + physiology*condition*clade, ExpDesign)
colnames(design)
##  [1] "physiologyFW"                           
##  [2] "physiologyM"                            
##  [3] "condition15_ppt"                        
##  [4] "cladeClade2"                            
##  [5] "cladeClade3"                            
##  [6] "physiologyM:condition15_ppt"            
##  [7] "physiologyM:cladeClade2"                
##  [8] "physiologyM:cladeClade3"                
##  [9] "condition15_ppt:cladeClade2"            
## [10] "condition15_ppt:cladeClade3"            
## [11] "physiologyM:condition15_ppt:cladeClade2"
## [12] "physiologyM:condition15_ppt:cladeClade3"
# check rank of matrix
#Matrix::rankMatrix( design )
#dim(design)

3 Filtering and Normalization

Genes with low expression across samples were dropped from the analysis using a conservative approach. The function filterByExpr was used on the raw counts matrix. For each condition_physiology group (regardless of species), each sample must have a minium count of 10, and a group minium total count of 100. This reduced the counts table to the following dimensions (genes x samples):

gene.names<-rownames(counts)
counts<-as.matrix(as.data.frame(sapply(counts, as.numeric)))
rownames(counts)<-gene.names
#class(counts)

keep<-filterByExpr(counts,design = design,group=condition_physiology,min.count = 10, min.total.count = 100)
counts.filt <- counts[keep,]
dim(counts.filt)
## [1] 21401    81
#write.csv(counts.filt,"../../../21k_counts_filt_30April2019.csv")

4 Genes of Interest

After filtration, we will check the counts matrix for the presence of several genes of interest. These genes have demonstrated responsiveness to salinity change from past studies.

gene Funhe Ensembl
zymogen granule membrane protein 16 Funhe2EKm029929 ENSFHEP00000007220.1
zymogen granule membrane protein 16 Funhe2EKm029931 ENSFHEP00000025841
solute carrier family 12 member 3-like (removed) Funhe2EKm006896 ENSFHEP00000009214
chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed) Funhe2EKm024148 ENSFHEP00000019510
ATP-sensitive inward rectifier potassium channel 1 Funhe2EKm001965 ENSFHEP00000015383
inward rectifier potassium channel 2 Funhe2EKm023780 ENSFHEP00000009753
aquaporin-3 ENSFHEP00000006725
cftr Funhe2EKm024501 ENSFHEP00000008393
polyamine-modulated factor 1-like Funhe2EKm031049 ENSFHEP00000013324
sodium/potassium/calcium exchanger 2 Funhe2EKm025497 ENSFHEP00000034177
septin-2B isoform X2 ENSFHEP00000015765
CLOCK-interacting pacemaker-like Funhe2EKm026846 ENSFHEP00000017303
vasopressin V2 receptor-like Funhe2EKm026721 ENSFHEP00000000036
sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1 Funhe2EKm001174 ENSFHEP00000031108
septin-2 Funhe2EKm012182 ENSFHEP00000016853
otopetrin-2 Funhe2EKm035427 ENSFHEP00000026411
claudin 8 Funhe2EKm037718 ENSFHEP00000006282
claudin 4 Funhe2EKm007149 ENSFHEP00000003908

If the Ensembl ID displays below, the gene is present in the whole data set and has not been filtered.

# ---------------------------
# Andrew Whitehead's genes of interest:
# ---------------------------

# zymogen granule membrane protein 16
# Funhe2EKm029929
# ENSFHEP00000007220.1
countsfilt <- counts.filt
countsfilt$row <- rownames(countsfilt)
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000007220.1"]
goi
## [1] "ENSFHEP00000007220.1"
# zymogen granule membrane protein 16
# Funhe2EKm029931
# ENSFHEP00000025841
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000025841"]
goi
## [1] "ENSFHEP00000025841"
# solute carrier family 12 member 3-like (removed) 
# Funhe2EKm006896
# ENSFHEP00000009214
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009214"]
goi
## [1] "ENSFHEP00000009214"
# chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed)
# Funhe2EKm024148
# ENSFHEP00000019510
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000019510"]
goi
## [1] "ENSFHEP00000019510"
# ATP-sensitive inward rectifier potassium channel 1 
# Funhe2EKm001965
# ENSFHEP00000015383
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015383"]
goi
## [1] "ENSFHEP00000015383"
# inward rectifier potassium channel 2
# Funhe2EKm023780
# ENSFHEP00000009753
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009753"]

# --------------------------------
# other salinity genes of interest
# --------------------------------
# ============================================
# aquaporin-3
# ENSFHEP00000006725
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006725"]
goi
## [1] "ENSFHEP00000006725"
# cftr
# Funhe2EKm024501
# ENSFHEP00000008393
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000008393"]
goi
## [1] "ENSFHEP00000008393"
# polyamine-modulated factor 1-like
# Funhe2EKm031049
# ENSFHEP00000013324
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000013324"]
goi
## [1] "ENSFHEP00000013324"
# sodium/potassium/calcium exchanger 2
# ENSFHEP00000034177
# Funhe2EKm025497
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000034177"]
goi
## character(0)
# septin-2B isoform X2
# ENSFHEP00000015765
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015765"]
goi
## [1] "ENSFHEP00000015765"
# CLOCK-interacting pacemaker-like
# ENSFHEP00000017303
# Funhe2EKm026846
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000017303"]
goi
## [1] "ENSFHEP00000017303"
# vasopressin V2 receptor-like
# Funhe2EKm026721
# ENSFHEP00000000036
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000000036"]
goi
## [1] "ENSFHEP00000000036"
# sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1
# ENSFHEP00000031108
# Funhe2EKm001174
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000031108"]
goi
## [1] "ENSFHEP00000031108"
# septin-2
# Funhe2EKm012182
# ENSFHEP00000016853
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000016853"]
goi
## [1] "ENSFHEP00000016853"
# otopetrin-2
# Funhe2EKm035427
# ENSFHEP00000026411
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000026411"]
goi
## [1] "ENSFHEP00000026411"
# claudin 8
# Funhe2EKm037718
# ENSFHEP00000006282
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006282"]
goi
## [1] "ENSFHEP00000006282"
# claudin 4
# ENSFHEP00000003908
# Funhe2EKm007149
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000003908"]
goi
## [1] "ENSFHEP00000003908"
all_goi<-c("ENSFHEP00000007220.1","ENSFHEP00000025841","ENSFHEP00000019510",
           "ENSFHEP00000015383","ENSFHEP00000009753","ENSFHEP00000006725","ENSFHEP00000008393",
           "ENSFHEP00000013324","ENSFHEP00000001609","ENSFHEP00000013324","ENSFHEP00000034177",
           "ENSFHEP00000015765","ENSFHEP00000017303","ENSFHEP00000000036","ENSFHEP00000031108",
           "ENSFHEP00000016853","ENSFHEP00000003908")

5 Exploratory Plots

Log counts before normalization:

# log counts before DE
boxplot(log(counts.filt+1), las = 2, main = "")

Log counts after cpm normalization

# ---------------

# DE analysis

# ---------------

genes = DGEList(count = counts.filt, group = condition_physiology_clade)
genes = calcNormFactors( genes )

# write normalized counts
dir <- "~/Documents/UCDavis/Whitehead/"
tmp <- as.data.frame(cpm(genes))
tmp$Ensembl <- rownames(tmp)
tmp <- dplyr::select(tmp, Ensembl, everything())
write.csv(tmp, file = file.path(dir, "normalized_counts.csv"), quote = F, row.names = F)

vobj = voom( genes, design, plot=TRUE)

lcpm <- cpm(genes$counts, log = TRUE)

# log counts after DE

boxplot(lcpm, las = 2, main = "")

plot(colSums(t(lcpm)))

Voom weights:

vwts <- voomWithQualityWeights(genes, design=design, normalization="quantile", plot=TRUE)

The random effects of species are taken into account with the duplicateCorrelation function, which estimates the correlation between species. This reflects the between-species variability. The higher the correlation (0-1.0), the higher the variability between species.

corfit <- duplicateCorrelation(vobj,design,block=ExpDesign$species)

corfit$consensus
## [1] 0.7596421
#[1] 0.751381

5.0.1 Individuals clustered by overall expression

counts_round<-round(counts.filt,digits=0)
dds <- DESeqDataSetFromMatrix(countData = counts_round,colData = ExpDesign,design = design)
## converting counts to integer mode
rld <- vst(dds, blind = FALSE,fitType='local')
sampleDists <- dist(t(assay(rld)))
df <- as.data.frame(colData(dds)[,c("physiology","condition","clade")])
sampleDistMatrix <- as.matrix( sampleDists )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors, annotation = df, show_rownames=F)

5.0.2 Individuals clustered by Top 100 genes

5.0.3 Individuals clustered by top 50 gene expression

5.0.4 PCA for overall expression

cowplot::plot_grid( plotPCA(rld, intgroup="condition"),
                    plotPCA(rld, intgroup="physiology"),
                    plotPCA(rld, intgroup="clade"),
                    plotPCA(rld, intgroup=c("clade","physiology","condition")),
                           align="c", ncol=2)

6 Contrasts

After running the lmFit function, which fits the linear model for each gene in the matrix and takes the random effects correlation into account, the resulting linear model fit is then used to compute the estimated coefficients and standard errors for a given set of contrasts:

fitRan <- lmFit(vobj,design,block=ExpDesign$species,correlation=corfit$consensus)
#colnames(coef(fitRan))
y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|clade*physiology ), reverse = T)
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$clade, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts

contrasts
## [1] "15_ppt_v_0.2_ppt_Clade1_FW" "15_ppt_v_0.2_ppt_Clade2_FW"
## [3] "15_ppt_v_0.2_ppt_Clade3_FW" "15_ppt_v_0.2_ppt_Clade1_M" 
## [5] "15_ppt_v_0.2_ppt_Clade2_M"  "15_ppt_v_0.2_ppt_Clade3_M"
tables <- lapply(contrasts, function(contr){
    #print(contr)
    cm <- contrast.matrix[contr,]
    ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
    cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
    tmp <- contrasts.fit(fitRan, contrasts = cm)
    tmp <- eBayes(tmp)
    tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
    #tmp3 <- tmp2
    #tmp3$row <- rownames(tmp3)
    #tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
    #tmp3 <- tmp3[order(tmp3$adj.P.Val),]
    filename <- paste0(contr, ".csv")
    #write.csv(tmp2, file = file.path(dir, filename),quote = F)
    tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
    header1 <- 6
    names(header1) <- paste0("Top 20 genes for ", contr)
    header2 <- 6
    names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
    header3 <- 6
    names(header3) <- paste0("Output file is ", filename)
    tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
    tab <- tab %>% kable_styling()
    return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}
)

7 Three-way contrasts

The number of genes significant for the three-way interaction of condition:physiology:clade:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")
Overview of results
Comparison Number of genes with adjusted P < 0.05
15_ppt_v_0.2_ppt_Clade1_FW 13
15_ppt_v_0.2_ppt_Clade2_FW 105
15_ppt_v_0.2_ppt_Clade3_FW 46
15_ppt_v_0.2_ppt_Clade1_M 533
15_ppt_v_0.2_ppt_Clade2_M 272
15_ppt_v_0.2_ppt_Clade3_M 105

8 Two-way contrasts

These genes demonstrate a conserved response to experimental condition across M or FW physiologies, regardless of clade.

y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|physiology ), reverse = T)
## NOTE: Results may be misleading due to involvement in interactions
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts

contrasts
## [1] "15_ppt_v_0.2_ppt_FW" "15_ppt_v_0.2_ppt_M"
tables <- lapply(contrasts, function(contr){
  #print(contr)
  cm <- contrast.matrix[contr,]
  ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
  cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
  tmp <- contrasts.fit(fitRan, contrasts = cm)
  tmp <- eBayes(tmp)
  tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
  #tmp3 <- tmp2
  #tmp3$row <- rownames(tmp3)
  #tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
  #tmp3 <- tmp3[order(tmp3$adj.P.Val),]
  filename <- paste0(contr, ".csv")
  #write.csv(tmp2, file = file.path(dir, filename),quote = F)
  tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
  header1 <- 6
  names(header1) <- paste0("Top 20 genes for ", contr)
  header2 <- 6
  names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
  header3 <- 6
  names(header3) <- paste0("Output file is ", filename)
  tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
  tab <- tab %>% kable_styling()
  return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}

)

The number of genes significant for the two-way interaction of condition:physiology, independent of clade:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")
Overview of results
Comparison Number of genes with adjusted P < 0.05
15_ppt_v_0.2_ppt_FW 71
15_ppt_v_0.2_ppt_M 616

9 Condition - main effect

These genes demonstrate a main effect of condition (15 ppt vs. 0.2 ppt), regardless of clade or physiology.

y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade, 
                data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition), reverse = T)
## NOTE: Results may be misleading due to involvement in interactions
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
rownames(contrast.matrix) <- contrasts

contrasts
## [1] "15_ppt_v_0.2_ppt"
tables <- lapply(contrasts, function(contr){
  #print(contr)
  cm <- contrast.matrix[contr,]
  ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
  cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
  tmp <- contrasts.fit(fitRan, contrasts = cm)
  tmp <- eBayes(tmp)
  tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
  #tmp3 <- tmp2
  #tmp3$row <- rownames(tmp3)
  #tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
  #tmp3 <- tmp3[order(tmp3$adj.P.Val),]
  filename <- paste0(contr, ".csv")
  #write.csv(tmp2, file = file.path(dir, filename),quote = F)
  tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
  header1 <- 6
  names(header1) <- paste0("Top 20 genes for ", contr)
  header2 <- 6
  names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
  header3 <- 6
  names(header3) <- paste0("Output file is ", filename)
  tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
  tab <- tab %>% kable_styling()
  return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}

)

The number of genes significant for the main effect condition:

sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
                     check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
  add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")
Overview of results
Comparison Number of genes with adjusted P < 0.05
15_ppt_v_0.2_ppt 207

10 Heatmaps

10.1 Marine-Clade 1 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.2 Freshwater-Clade 1 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.3 Marine-Clade 2 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.4 Freshwater-Clade 2 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.5 Marine-Clade 3 (three-way) response

rld <- log2(mean_norm_counts_ordered_M_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.6 Freshwater-Clade 3 (three-way) response

rld <- log2(mean_norm_counts_ordered_FW_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.6.1 Marine physiology (two-way) conserved response

rld <- log2(mean_norm_counts_ordered_M_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.6.2 Freshwater physiology (two-way) conserved response

rld <- log2(mean_norm_counts_ordered_FW_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

10.6.3 Condition response

rld <- log2(mean_norm_counts_ordered_condition_sig+1)
geneDists <- dist(mean_norm_counts_ordered_condition_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
         clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] data.table_1.12.2           emmeans_1.3.4              
##  [3] kableExtra_1.1.0            biomaRt_2.38.0             
##  [5] vsn_3.50.0                  lattice_0.20-38            
##  [7] gplots_3.0.1.1              edgeR_3.24.3               
##  [9] limma_3.38.3                DESeq2_1.22.2              
## [11] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
## [13] BiocParallel_1.16.6         matrixStats_0.54.0         
## [15] Biobase_2.42.0              GenomicRanges_1.34.0       
## [17] GenomeInfoDb_1.18.2         IRanges_2.16.0             
## [19] S4Vectors_0.20.1            BiocGenerics_0.28.0        
## [21] ggrepel_0.8.0               knitr_1.22                 
## [23] tidyr_0.8.3                 dplyr_0.8.0.1              
## [25] RColorBrewer_1.1-2          cowplot_0.9.4              
## [27] pheatmap_1.0.12             gtools_3.8.1               
## [29] ggplot2_3.1.1               MASS_7.3-51.4              
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       estimability_1.3       htmlTable_1.13.1      
##  [4] XVector_0.22.0         base64enc_0.1-3        rstudioapi_0.10       
##  [7] affyio_1.52.0          bit64_0.9-7            AnnotationDbi_1.44.0  
## [10] mvtnorm_1.0-10         xml2_1.2.0             splines_3.5.2         
## [13] geneplotter_1.60.0     Formula_1.2-3          jsonlite_1.6          
## [16] annotate_1.60.1        cluster_2.0.8          BiocManager_1.30.4    
## [19] readr_1.3.1            compiler_3.5.2         httr_1.4.0            
## [22] backports_1.1.4        assertthat_0.2.1       Matrix_1.2-17         
## [25] lazyeval_0.2.2         acepack_1.4.1          htmltools_0.3.6       
## [28] prettyunits_1.0.2      tools_3.5.2            coda_0.19-2           
## [31] gtable_0.3.0           glue_1.3.1             GenomeInfoDbData_1.2.0
## [34] affy_1.60.0            Rcpp_1.0.1             gdata_2.18.0          
## [37] preprocessCore_1.44.0  xfun_0.6               stringr_1.4.0         
## [40] rvest_0.3.2            statmod_1.4.30         XML_3.98-1.19         
## [43] zlibbioc_1.28.0        scales_1.0.0           hms_0.4.2             
## [46] yaml_2.2.0             memoise_1.1.0          gridExtra_2.3         
## [49] rpart_4.1-13           latticeExtra_0.6-28    stringi_1.4.3         
## [52] RSQLite_2.1.1          highr_0.8              genefilter_1.64.0     
## [55] checkmate_1.9.1        caTools_1.17.1.2       rlang_0.3.4           
## [58] pkgconfig_2.0.2        bitops_1.0-6           evaluate_0.13         
## [61] purrr_0.3.2            labeling_0.3           htmlwidgets_1.3       
## [64] bit_1.1-14             tidyselect_0.2.5       plyr_1.8.4            
## [67] magrittr_1.5           R6_2.4.0               Hmisc_4.2-0           
## [70] DBI_1.0.0              pillar_1.3.1           foreign_0.8-71        
## [73] withr_2.1.2            survival_2.44-1.1      RCurl_1.95-4.12       
## [76] nnet_7.3-12            tibble_2.1.1           crayon_1.3.4          
## [79] KernSmooth_2.23-15     rmarkdown_1.12         progress_1.2.0        
## [82] locfit_1.5-9.1         grid_3.5.2             blob_1.1.1            
## [85] digest_0.6.18          webshot_0.5.1          xtable_1.8-3          
## [88] munsell_0.5.0          viridisLite_0.3.0